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Abstract. Autonomous oscillations in biological systems may have a biochemical 
origin or result from an interplay between force-generating and visco-elastic elements. 
In molecular motor assemblies the force-generating elements are molecular engines 
and the visco-elastic elements are stiff cytoskeletal polymers. The physical mechanism 
leading to oscillations depends on the particular architecture of the assembly. Existing 
models can be grouped into two distinct categories: systems with a delayed force 
activation and anomalous force-velocity relations. We discuss these systems within 
phase plane analysis known from the theory of dynamic systems and by adopting 
methods from control theory, the Nyquist stability criterion. 
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1. Introduction 

Oscillations are ubiquitous in biological systems. Examples range from circadian 
rhythms associated with external periodicities to autonomous bio-chemical [HI] and 
mechano-chemical oscillators. Voltage oscillations in systems of ion channels embedded 
in a cell membrane are one of the most prominent examples for a biochemical oscillator. 
At the heart of its explanation within the Hodgkin-Huxley model there is a negative 
differential current-voltage relation for the ion channels. Such biochemical oscillators 
are by now well studied jSI]. Recently, there has been growing interest in mechano- 
chemical systems composed of visco-elastic biopolymer arrays and active molecular 
machines driven by the chemical cycle of ATP hydrolysis; these will be the focus of 
our contribution. We ask how the interplay between force-generation and energy-input 
by the molecular motors and the restoring forces and damping mechanisms provided by 
the biopolymer system may lead to oscillatory behavior. 

There are many examples of functional units in cellular systems whose main 
components are molecular machines and visco-elastic elements. The part list of auditory 
hair bundles contains stereocilia (elastic rods composed of bundles of stiff biopolymers F- 
actin), myosin (a force-generating molecular motor) and mechanosensitive ion-channels. 
Spontaneous oscillations in such systems are by now well documented experimentally [2Z] 
and explained within simple theoretical models fOlllSl- This oscillatory behavior of hair 
bundles provides a mechanism for active amplification of acoustic signals |7j. Another 
example are eucaryotic cilia and flagella which are used by many small organisms to 
swim. Here the main structural element are axonemes, which consist of a cylindrical 
arrangement of elastically linked microtubules (another example for a stiff biopolymer) 
and an assembly of dynein motors located between neighboring microtubules. Forces 
generated by the dynein motors induce relative sliding between microtubules which in 
turn results into bending of the axonemes . Oscillatory behavior in such systems is 
generic and intimately connected with a Hopf bifurcation jH] . Sarcomeres, the contractile 
units of muscle, are composed of two sets of protein filaments: myosin filaments 
(containing myosin motors) and actin filaments. The shortening of the sarcomere is 
achieved by the actin and myosin filaments sliding over one another. While the primary 
function of most types of muscle is to generate a unidirectional force, some muscle 
(notably the asynchronous, also known as fibrillar or myogenic insect flight muscle) 
are specialized on autonomous oscillatory contraction. The oscillation mechanism 
usually includes a number of regulatory proteins that constitute the delayed stretch- 
activation mechanism. However, oscillations were also observed in skeletal muscle 
fibrils, mostly under non-physiological conditions [S2I- Theoretical studies show that 
even skeletal muscle myosin could produce oscillatory motion ^1 ^1 which is 
normally not observable because the contributions of different sarcomeres cancel out. 
However, a sudden change in applied load can lead to synchronization of sarcomeres 
and the oscillations become observable for a certain period of time [15 . Assemblies of 
molecular motors and stiff biopolymers also play an essential role during cell division. 
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The duplicated genome is segregated into two daughter cells through the action of the 
mitotic spindle. It consists of fibres (microtubule bundles) radiating from two poles and 
meeting at the equator in the middle. The spindle poles have been found to oscillate 
jl6j . Recently it has been shown that these oscillations are again due to an interplay 
between force-generating and elastic elements [T7] . 

Here we will focus on a discussion of the generic features of oscillations in molecular 
motor assemblies to highlight the common physical mechanisms. For more system 
specific and detailed discussions we refer the reader to the literature. The defining 
elements of molecular motor assemblies are: (i) an external energy source driving 
the system towards a non-equilibrium steady state, (ii) viscoelastic elements providing 
restoring forces and damping mechanisms, (iii) possibly force- dependent biochemical 
reactions which allow for switching between different states of the molecular motor 
assembly. The particular interconnection between these components gives rise to the 
system response which may be controlled by the biochemical reaction rates. Defined as 
such there are obvious analogies to system and control theory which is a well established 
discipline in engineering ^21- In particular, linear system analysis is a most fruitful 
concepts which we will also employ to characterize the system's dynamic response. 

The mechanisms for oscillations in chemochemical models discussed in the literature 
can be grouped into two distinct categories. In delayed force activation systems an 
imposed displacement x is transduced into a force / with some delay time r (Debye 
relaxator). As such the system would evolve towards a time-independent steady state. 
Oscillatory response is achieved only after coupling the Debye relaxator to other dynamic 
elements with either a negative stiffness or some inertial (massive) load. Anomalous 
force-velocity relations are an alternative route towards oscillatory behavior. This means 
that an assembly of motors can move with two different velocities (e.g., one positive and 
one negative) under the same load If attached to an external spring, the motors 
will move forwards until the spring force exceeds a certain threshold, then they will 
switch to the other stable state and slide backwards, until another threshold is reached 
and so the cycle repeats 



2. Response Functions 

To assess the stability of an active mechanical system like an insect muscle attached to 
the wing or the dynein motor acting between two tubulin filaments in a flagellum, the 
response function can be defined the following way. We connect the active system to a 
mechanical actuator that keeps its length at the desired point of operation and at the 
same time imposes small oscillations with the amplitude xq around that point 

Ax(t) = Rexoe^'^* . (1) 

Fur sufficiently small oscillations the system responds with a force that oscillates with 
the same frequency 

(A/(t)) = Re/oe*-*. (2) 
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Here () denote an ensemble average over different realizations of the response of the 
stochastic system. Note that the sign in the exponent of the Fourier transform was 
chosen in the way that is most common in the literature on muscle, although it is 
opposite to the convention frequently used in physics. We define the force in the way 
that it is positive if the motor system is being pulled upon in the direction of positive 
X. In the limit of small amplitude Xq, the relation becomes linear 



with the frequency dependent response function or modulus G{uj). Here the real part 
of G corresponds to the elastic and the imaginary part to the dissipative response. For 
a spring with a spring constant k the response function is G{lj) = k and for a damping 
element ("dashpot") G{lj) = iwj. 

One may also encounter a situation where for a given force / one is observing the 
displacement x for various realizations of the system response, 



if the following two conditions are fulfilled. First, the system needs to be stable, so that 
the operational point is the same regardless whether the position or the force is imposed. 
Second, the systems are large enough such that variance of the stochastic variable is 
small compared to its mean. 

3. The Nyquist stability criterion 

The Nyquist stability criterion is a convenient tool to analyze the stability of a dynamical 
system if the complex response function G{uj) is known over the whole frequency range. 
It states that the system is dynamically unstable if the plot of lmG{uj) vs. ReG(ci;) 
(Nyquist plot, sometimes also called Cole-Cole plot) encircles the coordinate origin in 
clockwise direction. It is frequently used in engineering, especially for the stability 
analysis of closed feedback loops in electrical circuits [T2j . 

To derive it we first note that the system is dynamically unstable if it starts 
oscillating with a growing amplitude in the absence of a force, (/) = 0. This is the 
case if there exists an eigenmode, i.e., a complex solution of the equation {f{xo, u)) = 
with Imcj < 0, since then the corresponding solution Ax{t) oc e*"^* obviously diverges for 
t +00. This means that G{uj) has a zero and x(^) has a pole in the lower half-plane 
with Imuj < 0. The Nyquist criterion is based on the fact that G{uj) is an analytical 
complex function and utilizes Cauchy's theorem to relate the number of zeros in the 
negative imaginary uj half-plane to the contour of G{u!), evaluated for real values of u. 

According to Cauchy's theorem, the change of phase along a closed contour in 
clockwise direction in the complex u plane equals 



/o = G{uj)xo , 



(3) 




(4) 



The two response functions become equivalent 




(5) 
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where is the number of encircled zeros of G{uj), Np the number of poles and Uz^k 
and rip their multiplicities. For example, the function G{uj) = iwy has a single zero 
{Nz = 1) with multiplicity 1 (n^^i = 1) at u; = 0. Accordingly, it changes its phase by 
— 27r on a clockwise circle around the origin. The function G{uj) = —muj'^, on the other 
hand, has a zero with multiplicity n^^i = 2 at u = and it changes its phase by —4n on 
the same circle. 

We apply Cauchy's theorem to a path that runs along the real u axis from uj = —oo 
to uj = +00 and close it in the negative imaginary plane as shown in Fig.^. We further 
note that in a real system G{uj) can have no poles in the negative imaginary half-plane, 
because they would imply an infinite force response to a finite displacement growing in 
time. Therefore, any phase change along that path indicates a zero with Imuj < 0, and 
therefore a dynamical instability. 
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Figure 1. a) The integration path for Eq. Elin the complex cu plane. Encircled zeros 
(marked "x") correspond to unstable eigenmodes. b) G in the complex plane, with 
as a parameter. The points denoted as A,B and C correspond to a negative, and 
a positive value of u. The curve shown has a typical shape measured on insect flight 
muscle SSj. 



Now we can look at the same path in the complex G plane, i.e., lmG{uj) plotted 
against ReG(u;), where u acts as a curve parameter. The strand along the real axis 
in the u plane corresponds to the (ReG{uj),lmG{u!)) curve, with u taking real values 
values from — oo to oo. In case G{—oo) ^ G'(oo), we close the path in the positive 
real half-plane (Fig. ^), which corresponds to the closure in the Imcj < half-plane. 
A change of phase of G along the path in the uj plane by a multiple of —In directly 
represents the same number of origin encirclements in the G plane. Therefore, it follows 
from (jH)) that the complex G{uS) curve encircles the coordinate origin Yl!\'' '^z,k times in 
clockwise direction. 

Note that the curve is symmetric with respect to the real axis. This becomes 
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evident from the following consideration. Because G{t) is a real function, its Fourier 
transformed has to fulfill the following symmetry relation: 



For real u values, this means G*{uj) = G{—uj), and for each G value its complex 
conjugate is part of the curve as well. 

An additional consequence of the Nyquist criterion is that the origin can only be 
encircled in clockwise direction, because G contains zeros but no poles in the relevant 
region. 

4. Delayed force activation 

One type of models, often used to explain oscillations generated by molecular motors 
involves a delayed force activation mechanism. This means that when the system is 
displaced by the distance x in one direction, it develops a restoring force /, however not 
instantaneously, but according to a differential equation like: 



Here r represents the time constant with which the force reacts to a position change 
and Ka an effective steady state stiffness of the active system. There are many different 
physical mechanisms that can produce a delayed force response (Fig. |21): 

(i) Stretch-activation of insect flight muscle: a number of regulatory proteins 
facilitate the activation of myosin motors following a mechanical stretch (a relative 
displacement of thin filaments relative to thick) [HHl ITT| IT]. After the activation, 
the development of force does not follow instantaneously, but with a certain time 
constant, typically with the attachment rate of motors. 

(ii) In a sarcomere a position change can bring myosin heads closer to the accessible 
actin sites jSHEH]- This mechanism has been proposed for insect flight muscle jH], 
although experiments have later refuted it in many cases |86] . 

(iii) The switching of bound myosin heads between two conformations with different 
detachment rates can lead to an effective stretch-dependent activation IH]. 

(iv) The opening of a mechanically gated ion channel can cause the influx of calcium 
ions, which then promote channel reclosure. This is called the fast adaptation 
mechanism in auditory hair cells and can generate spontaneous hair bundle 
oscillations jHH 1201 GDI 1121 • ^^^^ ^^^^ verge of spontaneous oscillations 
(i.e., close to a Hopf bifurcation) it has the best sensitivity, dynamic range and 
frequency selectivity [Z]. 

(v) Calcium ions, entering the ion channel, can reduce the force generated by adaptation 
motors of type myosin IC P], a mechanism known as slow adaptation in auditory 
hair cells [2ll|22l|2Hl. 




(7) 



df 



{Kax - f)/T . 



(8) 



dt 



Oscillations in molecular motor assemblies 



7 




Figure 2. Different biological mechanism that lead to delayed force activation: i) 
Stretch activation in insect flight muscle, after [313 . The applied stretch activates 
myosin mechanically through conformational changes involving flightin, GST-2, 
troponin-H and tropomyosin, ii) The alternative match-mismatch mechanism, where 
the myosin heads are activated when they are brought to appropriate positions along 
the actin helix, iii) The mechanism involving switching of myosin heads between two 
bound conformations. A stretch brings more myosin heads into the pre-powerstroke 
state, which is less likely to detach from actin, and therefore the total number of 
active heads increases [101 . iv) The fast adaptation mechanism in auditory hair cells. 
When the transduction channel opens, calcium ions enter and promote the reclosure 
of the channel, v) The slow adaptation mechanism in auditory hair cells: calcium ions 
entering the cells reduce the force generated by myosin motors, therefore they slip back 
and the channel is more likely to close again. 



(vi) Some models for flagellar dynein are based on a strain-controlled [2^] or curvature- 
controlled [H] activation. 

The solution of (jHl) can be expressed with a Green's function G'yi(t), so that 

f(t) = f GA{t - t')x{t')dt' with GA{t) = — e'*/" . (9) 

Such behavior is well known from many branches in physics under the term Debye 
relaxator. The corresponding response function 

Ga{uj) = (10) 

1 + lUT 
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is represented by a circle in a Nyquist plot (see Fig. |21 (a)). It does have a negative 
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Figure 3. Nyquist plots (also called Cole-Cole plots or vector modulus plots) of 
the response function of a system with delayed force activation, a) The response 
function of the delayed force activation alone, given by Eq. H1U|) . b) Response 
function of the delayed force activation, combined with a damping element, Ga + Gd 
(spring K — 5Ka, dashpot with a constant j/K = 100/r). c) System as in (b), 
additionally combined with an element with negative effective elasticity; Ga + Gd + Gn 
[Ki = --Q.2Ka)- d) Response function Ga + Gi of the system with delayed force 
activation, combined with a damped inertial oscillator (wq = 10/r, K = 2Ka, 
/3=10/r). 

imaginary part, however, the real part is always positive. Hence the Nyquist plot of 
Ga{uj) cannot encircle the coordinate origin and thereby satisfy the Nyquist criterion for 
an instability. The stationary state of the Debye relaxator simply is a time-independent 
displacement x* = f/Kj^. Note, that a dynamics with a one-dimensional phase space, 
i.e. first order differential equations of the form df/dt = x{f), is either monotonic or 
constant Only if the dynamic variable is on a circle or if the phase space has more 
than two dimensions oscillations emerge as possible solutions. 
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One way to make the system unstable is to connect the active mechanism with an 
element that shifts the Nyquist plot towards negative values on the real axis; for an 
illustration of the type of connection see Fig. IHl This can be trivially achieved through 
an element with a negative effective stiffness, or, as we will show below, through an 
inertial load. 

Figure El shows the response function of the delayed activation system (fTU)) . 
combined with different elements. The first is an element consisting of a spring (constant 
K) and a dashpot (damping coefficient 7) in series. Their combined response function 
reads 



A second element we connect to the motor assembly is an element with negative effective 
elasticity. 



The third element is a damped inertial oscillator (a spring connected to a mass, with an 
additional damping element connecting both to a fixed point). The response function 
of such an element reads 



with Uq = K/m and (3 = 7/m. 

Negative stiffness can be generated by the gating compliance of ion channels [^1^ 
or by motors switching between two states jTSl EH] • Inertial load is the likely mechanism 
of ensuring an instability in insect flight muscle [H] . Also in many types of hair cells, 
no negative stiffness has been observed, but as the hair bundle is often attached to the 
tectorial membrane, inertial load could play an important role as well 

5. Anomalous force-velocity relations 

A class of models that were frequently discussed in relation to dynamical instabilities in 
systems of molecular motors is based on an anomalous force- velocity relation, meaning 
that a group of coupled motors produces a force that displays a region with negative 
slope as a function of velocity. Models for dynein 0, myosin propulsion by actin 
polymerization |3j , mitotic spindle oscillations [T7] as well as generic motors [23 121] 
were based on this mechanism. 

A common feature of such models is translational invariance, meaning that the 
force fit) does not depend on the position x, but only on the velocity v and the history 
thereof. The translational invariance can be achieved even on periodic tracks (like 
actin filaments) if the motors (myosin heads) are arranged incommensurate to the track 
periodicity [SB]- 

An instructive example of a system displaying an anomalous force- velocity relation 
is a two-state cross-bridge model for motors like myosin. In this model, motors can 






(12) 




(13) 
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attach to actin in a forward-leaning position with rate r^, then quickly undergo a 
conformational change, which makes them strained. When the filaments moves, this 
strain changes with time (for forward-running motors, it decreases), and when the motor 
eventually detaches with a strain-dependent rate r^{^) the cycle can repeat. Each motor 
generates a force k^, proportional to its strain ^. A full description of this system requires 
Master equations for the probability of a motor for being in the detached state ($d) 
and for the probability density of a motor for being in the attached state with strain 

{dt - vd^) <l>,(e, t) = $d(t)raP(0 - '^.i^, t)rM) (14a) 

/oo 
de$a(e,t)rd(o, im 

Both probabilities are not independent but normalized such that $d + $a(0'^^ = ^ 
gives the total number of motors. P{C,), normalized as P{^)dC, = 1 is the distribution 
of strains on newly attached motors, with an expectation value d = /f^-P(OCc?^- A 
general analytical solution for the stationary state of these equations is given in • 

To illustrate the key physical features of this model we neglect the strain- 
dependence of the attachment probability and assume that all attached motors have 
the same strain y(t), 

M^,t) = n{t)5{^-yit)). (15) 

Then the number of attached motors n{t) is related to the number of detached motors 
by ^d{t) = N — n{t). Upon inserting these relations into Eqs ll4a^ b , multiplying Eq ll4flJ 
with ^ and integrating over all values of ^ one obtains a simplified set of dynamic 
equations for y(t) and n(t): 

N — n 

y = [d-y] r^-v (16a) 

n 

h = [N — n]ra — nrA{y) (166) 

The external force on the system, which equals the total force produced by the motors, 
is given as 

/ = -nky (17) 

where k denotes the spring constant of each motor. These are nonlinear flow equations 
with a two-dimensional phase space consisting of the strain y{t) and the number density 
n{t) of the attached motors. For a constant velocity f , they always have a stable 
stationary solution. However, if we take the system consisting of the group of motors, 
coupled to an external spring with a spring constant so that the force equilibrium 
states / + Kx = 0, we obtain a new system of two nonlinear equations for two 
independent variables (e.g., x and n). Depending on the parameter values, they can 
have a stable fixed point or a limit cycle (oscillatory steady state). These two examples 
are shown in Fig. |31 

In the following we will again use the stability analysis of linear response functions 
to characterize the system. In the absence of an external velocity {v = 0), Eqs. 116 a[ b 
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Figure 4. Numerical solution of the equations of motion p6al b') for the system of 
motors, coupled to a spring with a constant K = 0.4 x Nk (a) and K = 0.3 x Nk (b). 
The force-dependent detachment rate was chosen as rd{y) — Va exp(2(x/d)^ — 2). The 
system has a stable fixed point in the first case and an unstable one, surrounded by a 
limit cycle, in the second case. 



have a stable stationary solution with i/q = d and no 
we can linearize the flow equations to 

y = - rd{d)Ay - v 

h = - nor'^{d)Ay - {r^ + raid)) An 



r-a+rd{rf) 



A^. For small values of v 

(18a) 
(186) 



where Ay = y — y^ and An = n — n^. These equations have two negative eigenvalues, 
— Ai = —rd{d) and — A2 = — Ta — rd((i). The solution of the coupled system then has the 
form 

Ay{t) = - r dt'e^^'^'-''^v{t') (19a) 



An{t) 



nor'M) 



dt' 



'X-2{t-t') 



v{t') 



The force can be determined as 



f{t) = -nky = fo + 



dt'knn 



(196) 



v{t') . (20) 



Where a 



is a dimensionless coefficient. With a more precise solution, or with 



more complex model equations, the solution could contain more than two relaxation 
constants, but its basic form would remain. Inserting a stationary velocity {v{t) = const) 
in Eq. pUjl one can see that the system displays an anomalous stationary force velocity 
relation if a is sufficiently large that 

" " (21) 



Ai A2 

The response function in Fourier space follows from Eq. ()20|) : 



G{lu) = no/c(l — a) 



Ai + iuj 



+ n^ka 



A2 + ioj 



(22) 
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In the limit — oo, the response is always that of the elastic elements involved, 
G = n^k. Figure El shows the resulting Nyquist plot for two scenarios: with a normal 
and an anomalous force- velocity relation. In case of an anomalous force velocity relation 
the curve encircles the origin and the system is dynamically unstable. 




Figure 5. Nyquist plots (Cole-Cole plots, vector modulus plots) of the response 
function G{Ld) H22|) for a = (a) and a = 5 (b). Other parameter values: rd(d) = ra, 
implying A2 = 2Ai. Negative uj values are shown with dashed lines. The system in (a) 
is dynamically stable, because the curve does not encircle the origin, while that in (b) 
oscillates spontaneously. 



If the assembly of motors is coupled to an elastic element (spring constant K) , the 
curve in the Nyquist plot is shifted to the right by the amount K. There is a critical 
value of K at which the oscillations stop. In the example shown in Fig. |3J this value 
would lie between those used in both diagrams. At the transition point the system 
exhibits a Hopf bifurcation pij I17j. 

In summary, based on the foregoing discussion, we feel that a combination of 
methods from the engineering sciences like control theory and methods from the theory 
of dynamic systems may be fruitful for future analysis of molecular motor assemblies or 
other functional units in cellular systems. 
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